Chapter 6: Model critique and evidence consistency in random effects meta-analysis
Bob
2025-02-28
library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots') # caterplot
library('dplyr')criticism_results
Random effects model
# Model
model_code <- "model
{
for(i in 1:32){
r[i] ~ dbin(p[i], n[i]) #Likelihood
#Log-odds for control (t[i]=1) and treatment (t[i]=2) groups
logit(p[i]) <- mu[s[i]] + delta[s[i]]*equals(t[i],2)
#Deviance contribution
rhat[i] <- p[i] * n[i] # expected value of the numerators
dev[i] <- 2 * (r[i] * (log(r[i])-log(rhat[i])) + (n[i]-r[i]) * (log(n[i]-r[i]) - log(n[i]-rhat[i])))
}
for (j in 1: 16){
#Priors for baseline effects
mu[j] ~ dnorm(0,1.0E-6)
#Hierarchical random effects model for treatment effects
delta[j] ~ dnorm(d, prec)
}
or<-exp(d) #Population odds ratio
d ~ dnorm(0.0,1.0E-6) #Prior for population treatment effect
tau ~ dnorm(0.0,1.0E-6)I(0,) #Prior for between studies sd
prec<- 1/(tau*tau)
delta.new~dnorm(d,prec) #Replicate log OR for prediction
# For plotting purposes put population treatment effect & predicted effect from new study
# in elements 19 & 20 respectively
delta[19] <- d
delta[20] <- delta.new
#Total Deviance
resdev <- sum(dev[])
}
" %>% textConnection
# Data
data <- "s t r n
1 1 2 36
1 2 1 40
2 1 23 135
2 2 9 135
3 1 7 200
3 2 2 200
4 1 1 46
4 2 1 48
5 1 8 148
5 2 10 150
6 1 9 56
6 2 1 59
7 1 3 23
7 2 1 25
8 1 1 21
8 2 0 22
9 1 11 75
9 2 6 76
10 1 7 27
10 2 1 27
11 1 12 80
11 2 2 89
12 1 13 33
12 2 5 23
13 1 8 122
13 2 4 130
14 1 118 1157
14 2 90 1159
15 1 17 108
15 2 4 107
16 1 2103 29039
16 2 2216 29011" %>% read.delim(text=.) %>% as.list
# Starting/initial values
initial_values <- list(
list(d = 0, tau=1,mu = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
delta = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,NA,NA,NA,NA),
delta.new=0
),
list(d = -1, tau=2, mu = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
delta = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,NA,NA,NA,NA),
delta.new=0
),
list(d = 1, tau=0.5,mu = c(-1,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1),
delta = c(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,NA, NA,NA,NA),
delta.new=0
)
)
parameters_to_save <- c("d","delta", "dev", "or", "prec", "resdev", "tau")
results <- jags(data = data,
inits = initial_values,
parameters.to.save = parameters_to_save,
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 32
## Unobserved stochastic nodes: 35
## Total graph size: 668
##
## Initializing model
##
## Adaptive phase, 100 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 3 chains
##
##
## Sampling from joint posterior, 30000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.summary(results)## mean sd 2.5% 25% 50%
## d -0.90404672 0.27436264 -1.498971e+00 -1.07040605 -0.88628981
## delta[1] -0.89421922 0.66465245 -2.287918e+00 -1.29759367 -0.86638632
## delta[2] -1.02165360 0.36495353 -1.774609e+00 -1.26019716 -1.00735616
## delta[3] -1.09008801 0.56090320 -2.294945e+00 -1.44025620 -1.05673516
## delta[4] -0.73808237 0.69647378 -2.141689e+00 -1.16567311 -0.73457252
## delta[5] -0.16606609 0.42365551 -9.751880e-01 -0.45184314 -0.17458763
## delta[6] -1.49490200 0.62043872 -2.861497e+00 -1.86512747 -1.44036861
## delta[7] -1.02280997 0.65378130 -2.418614e+00 -1.42283670 -0.98497706
## delta[8] -1.07187369 0.77318832 -2.769949e+00 -1.51865204 -1.01682417
## delta[9] -0.77681697 0.43550751 -1.659283e+00 -1.06202548 -0.76817809
## delta[10] -1.37371519 0.63799404 -2.791584e+00 -1.75007628 -1.31932819
## delta[11] -1.47272171 0.55423410 -2.671270e+00 -1.81654983 -1.43242253
## delta[12] -0.88374776 0.47939189 -1.867594e+00 -1.19062217 -0.86933493
## delta[13] -0.85593265 0.48421207 -1.854932e+00 -1.16352476 -0.84337894
## delta[14] -0.33014205 0.14288969 -6.113890e-01 -0.42544151 -0.32966958
## delta[15] -1.32469315 0.46513230 -2.321688e+00 -1.61673391 -1.29823410
## delta[16] 0.05515721 0.03156944 -7.050354e-03 0.03395599 0.05534603
## delta[19] -0.90404672 0.27436264 -1.498971e+00 -1.07040605 -0.88628981
## delta[20] -0.90453094 0.81484732 -2.603124e+00 -1.37482915 -0.87563495
## dev[1] 0.78735352 1.10807900 7.468822e-04 0.08096559 0.36077211
## dev[2] 0.53477621 0.74678142 5.076015e-04 0.05533416 0.24874425
## dev[3] 0.94435397 1.35321453 8.452687e-04 0.09403138 0.42569557
## dev[4] 0.83235452 1.17899333 8.535195e-04 0.08569379 0.37932149
## dev[5] 0.95839621 1.32632842 9.722126e-04 0.09718985 0.44133241
## dev[6] 0.62435730 0.90658331 5.857951e-04 0.06104878 0.27529279
## dev[7] 0.67774845 0.97127379 6.724926e-04 0.06838231 0.30640112
## dev[8] 0.83807158 1.05226980 1.056953e-03 0.10464266 0.44580652
## dev[9] 1.05947697 1.46884244 1.032915e-03 0.10908266 0.49171403
## dev[10] 1.32284465 1.72753609 1.432552e-03 0.15064627 0.66926838
## dev[11] 1.25017494 1.68448128 1.276456e-03 0.13541751 0.60026459
## dev[12] 1.24572747 1.49305217 1.828778e-03 0.18141763 0.71496885
## dev[13] 0.88625192 1.23092127 9.562462e-04 0.09319689 0.40702141
## dev[14] 0.50900376 0.74751540 5.436790e-04 0.05102252 0.22750336
## dev[15] 1.30184239 1.72533552 1.370873e-03 0.14451621 0.63558217
## dev[16] 0.59705941 0.72151153 1.051376e-02 0.13253888 0.34963979
## dev[17] 0.86366925 1.22443981 8.345482e-04 0.08905054 0.39617209
## dev[18] 0.80133365 1.12969110 7.677557e-04 0.08175295 0.36657847
## dev[19] 1.17542425 1.59741712 1.165849e-03 0.12716474 0.55844176
## dev[20] 1.00506404 1.29062978 1.283094e-03 0.12342378 0.52457885
## dev[21] 1.17367299 1.60919120 1.208699e-03 0.12184268 0.54465025
## dev[22] 1.08040497 1.41836296 1.304302e-03 0.12302133 0.53539430
## dev[23] 0.87511817 1.23982980 8.825754e-04 0.08751937 0.39767095
## dev[24] 0.72395802 1.03681573 7.501261e-04 0.07405602 0.32310103
## dev[25] 0.86180987 1.20707065 8.809494e-04 0.08892404 0.39873539
## dev[26] 0.73426610 1.03158775 7.017103e-04 0.07463550 0.33391356
## dev[27] 0.98315903 1.37606510 9.880295e-04 0.09814564 0.44760834
## dev[28] 1.00101638 1.40784023 1.045017e-03 0.10143130 0.46040310
## dev[29] 1.05414925 1.47132647 1.014276e-03 0.10762229 0.48830996
## dev[30] 0.86773712 1.21393355 8.386151e-04 0.08883158 0.40525211
## dev[31] 0.99740841 1.40256148 8.443172e-04 0.09985963 0.45649128
## dev[32] 1.00157199 1.43044302 9.458723e-04 0.10067012 0.45143105
## or 0.41994720 0.11237226 2.233598e-01 0.34286927 0.41218219
## prec 2.69476693 2.08145531 5.770277e-01 1.36824353 2.13544447
## resdev 29.56955673 7.36229928 1.691502e+01 24.35465133 28.97018614
## tau 0.72180634 0.24855531 3.507458e-01 0.54656991 0.68431468
## deviance 147.35066741 7.36229928 1.346961e+02 142.13576201 146.75129682
## 75% 97.5% Rhat n.eff overlap0 f
## d -0.72026194 -0.41369940 1.000367 6562 0 0.9996000
## delta[1] -0.46647224 0.36237499 1.000168 22433 1 0.9224444
## delta[2] -0.77127144 -0.33524492 1.000059 27851 0 0.9984222
## delta[3] -0.70827617 -0.06770754 1.000092 23757 0 0.9817111
## delta[4] -0.30514920 0.66051211 1.000119 28083 1 0.8702444
## delta[5] 0.10872347 0.69350707 1.000096 21499 1 0.6610667
## delta[6] -1.06140024 -0.43354681 1.000136 13541 0 0.9978222
## delta[7] -0.58746164 0.18071606 1.000092 20868 1 0.9530889
## delta[8] -0.57100878 0.33369052 1.000516 4712 1 0.9351556
## delta[9] -0.48654131 0.06232654 1.000087 24886 1 0.9651778
## delta[10] -0.93465273 -0.27016972 1.000175 15661 0 0.9931778
## delta[11] -1.08517635 -0.50605425 1.000231 8485 0 0.9993556
## delta[12] -0.56537799 0.03262513 1.000099 18559 1 0.9706889
## delta[13] -0.53276981 0.06231826 1.000087 22341 1 0.9656667
## delta[14] -0.23406314 -0.04927802 1.000049 27370 0 0.9897111
## delta[15] -1.00343298 -0.48242332 1.000189 11814 0 0.9994000
## delta[16] 0.07633159 0.11704247 1.000173 16750 1 0.9588222
## delta[19] -0.72026194 -0.41369940 1.000367 6562 0 0.9996000
## delta[20] -0.41118986 0.67863958 1.000122 13479 1 0.8905778
## dev[1] 1.04026749 3.96167704 1.000611 17062 0 1.0000000
## dev[2] 0.71740205 2.61874164 1.000102 45000 0 1.0000000
## dev[3] 1.23792065 4.77623317 1.000071 45000 0 1.0000000
## dev[4] 1.10038789 4.19497627 1.000507 13008 0 1.0000000
## dev[5] 1.28986786 4.74213813 1.000547 7444 0 1.0000000
## dev[6] 0.81469641 3.17465660 1.000175 45000 0 1.0000000
## dev[7] 0.89285030 3.40195592 1.000519 45000 0 1.0000000
## dev[8] 1.18375576 3.77322027 1.000084 45000 0 1.0000000
## dev[9] 1.41699205 5.22151535 1.000591 13225 0 1.0000000
## dev[10] 1.84536534 6.16032135 1.000074 45000 0 1.0000000
## dev[11] 1.69004286 5.99292284 1.000188 12600 0 1.0000000
## dev[12] 1.77348433 5.38355679 1.000277 15380 0 1.0000000
## dev[13] 1.17966811 4.43787036 1.000463 29021 0 1.0000000
## dev[14] 0.65633988 2.62732306 1.000064 45000 0 1.0000000
## dev[15] 1.78535184 6.16953446 1.000023 45000 0 1.0000000
## dev[16] 0.78553457 2.60474052 1.000728 11160 0 1.0000000
## dev[17] 1.12871020 4.33454997 1.000350 12311 0 1.0000000
## dev[18] 1.06478149 3.97434027 1.000044 45000 0 1.0000000
## dev[19] 1.58625645 5.70697982 1.000181 35523 0 1.0000000
## dev[20] 1.40510905 4.58963068 1.000010 45000 0 1.0000000
## dev[21] 1.59432758 5.78060127 1.000100 45000 0 1.0000000
## dev[22] 1.49878815 5.02284087 1.000160 17824 0 1.0000000
## dev[23] 1.16384069 4.39252972 1.000457 16928 0 1.0000000
## dev[24] 0.94399265 3.71613117 1.000070 45000 0 1.0000000
## dev[25] 1.14929951 4.23245015 1.000613 17251 0 1.0000000
## dev[26] 0.98274602 3.69896492 1.000005 45000 0 1.0000000
## dev[27] 1.31463849 4.94617287 1.000023 45000 0 1.0000000
## dev[28] 1.31390095 5.03227690 1.000184 24879 0 1.0000000
## dev[29] 1.39679081 5.25197603 1.000063 45000 0 1.0000000
## dev[30] 1.15315016 4.28566326 1.000871 11616 0 1.0000000
## dev[31] 1.33160686 4.95317272 1.000385 29059 0 1.0000000
## dev[32] 1.32489573 5.12587001 1.000107 45000 0 1.0000000
## or 0.48662477 0.66119968 1.000227 8105 0 1.0000000
## prec 3.34740738 8.12858564 1.000797 8125 0 1.0000000
## resdev 34.13914734 45.68301988 1.000171 11073 0 1.0000000
## tau 0.85490587 1.31644180 1.000415 7104 0 1.0000000
## deviance 151.92025802 163.46413056 1.000171 11073 0 1.0000000plot(results)
